#!D:/Program Files/Python37/python
# -*- coding: utf-8 -*-
# @Time : 2020/7/9 10:55
# @Author : Xingzhi Zhao
# @File : core.py
# @Software : PyCharm

import math


def leastSquareLinearFit (x: list, y: list, num: int):
	"""最小二乘法算法实现"""
	sum_x2: float = 0.0
	sum_y2: float = 0.0
	sum_y: float = 0.0
	sum_x: float = 0.0
	sum_xy: float = 0.0

	for i in range(num):
		sum_x2 += x[i] * x[i]
		sum_y2 += y[i] * y[i]

		sum_x += x[i]
		sum_y += y[i]

		sum_xy += x[i] * y[i]

	ave_x = sum_x / num
	ave_y = sum_y / num

	a = (num * sum_xy - sum_x * sum_y) / (num * sum_x2 - sum_x * sum_x)
	b = (sum_x2 * sum_y - sum_x * sum_xy) / (num * sum_x2 - sum_x * sum_x)

	r = (sum_xy - num * (ave_x) * (ave_y)) / math.sqrt(
		(sum_x2 - num * (ave_x)**2) * (sum_y2 - num * (ave_y)**2))
	r2 = r * r

	return a, b, r2


def maxR2 (x: list, y: list, num: int):
	"""
	计算当r^2最大时b的值
	:param x: x轴自变量
	:type x: list
	:param y: y轴因变量
	:type y: list
	:param num: 列表长度
	:type num: int
	:return: α，β，最大的R^2及其对应的b
	:rtype: int
	"""
	b = 0.80
	max_r2 = 0
	lg = []
	t_bAdd1 = []
	alpha_fix, beta_fix, r2, b_fix = None, None, None, None
	while b <= 1.30:
		# print('b=', b)
		for i in range(num):
			lg.append(math.log10(y[i] / (x[i]**b)))
			t_bAdd1.append(x[i]**(b + 1))

		# print(lg, '\s', t_bAdd1)
		alpha, beta, r2 = leastSquareLinearFit(t_bAdd1, lg, len(x))
		# print('a=', alpha, 'b=', beta, 'r2=', r2)
		if max_r2 <= r2:
			max_r2 = r2
			alpha_fix = alpha
			beta_fix = beta
			b_fix = b
		b += 0.01
		lg = []
		t_bAdd1 = []
	print("alpha=", alpha_fix, "beta=", beta_fix, "max_r2=", max_r2, "b_fix=", b_fix)
	return alpha_fix, beta_fix, max_r2, b_fix


def getAandB (np_list: list, Lp_list: list, num: int):
	A, B, r2 = 0.0, 0.0, 0.0
	if len(np_list) == len(Lp_list):
		logLp_list = []
		for lp in Lp_list:
			logLp_list.append(math.log10(lp))
		B, A, r2 = leastSquareLinearFit(np_list, logLp_list, num)
	return A, B, r2


def getQoTest (t_list: list, num, a, b, c):
	QoTest_list = []
	for i in range(num):
		QoTest_list.append(a * t_list[i] * math.exp(-(t_list[i]**(b + 1) / c)))
	return QoTest_list


def getNpTest (t_list: list, num, a, b, c):
	NpTest_list = []
	for i in range(num):
		NpTest_list.append(a * c / (b + 1) * (1 - math.exp(-(t_list[i]**(b + 1)) / c)))
	return NpTest_list


def getFwTest (NpTest_list: list, num, A, B):
	FwTest_list = []
	for i in range(num):
		FwTest_list.append(1 - math.pow(10, -(A + B * NpTest_list[i] + math.log10(2.303 * B))))
	return FwTest_list


def getQmaxAndtm (a, b, c):
	Qmax = a * ((b * c / (b + 1))**(b / (b + 1))) * math.exp(-(b / b + 1))
	tm = b * c / (b + 1)**(1 / (b + 1))
	return Qmax, tm


def DataProcessing (t, Qo, Np, Lp, t_start=1, Lp_start=6):
	"""参数计算函数"""
	x = t[t_start:]
	y = Qo[t_start:]
	LpSample = Lp[Lp_start:]
	NpSample = Np[Lp_start:]

	alpha, beta, r2, b_fix = None, None, None, None
	if len(x) == len(y):
		beta, alpha, r2, b_fix = maxR2(x, y, len(x))
	else:
		print('length of x and y isn\'t equal.')
	print('alpha=', alpha, 'beta=', beta, 'r2=', r2, 'b=', b_fix)
	a = 10**alpha
	c = -1 / (2.303 * beta)
	Nr = a * c / (b_fix + 1)
	A, B, ABR2 = getAandB(NpSample, LpSample, len(NpSample))
	Qmax, tm = getQmaxAndtm(a, b_fix, c)
	Fw_limit = 1 - 10**(-(A + B * Nr + math.log10(2.303 * B)))
	return 0, A, B, ABR2, alpha, beta, r2, a, b_fix, c, Fw_limit, Qmax, tm, Nr


def predictionParameters (t_pre: int, a, b_fix, c, A, B):
	"""编写报告用的预测函数"""
	t = []
	for i in range(t_pre):
		t.append(i + 1)
	QoTest = getQoTest(t, len(t), a, b_fix, c)
	NpTest = getNpTest(t, len(t), a, b_fix, c)
	FwTest = getFwTest(NpTest, len(NpTest), A, B)

	print('\nQoTest =', QoTest, '\nNpTest =', NpTest, '\nFwTest =', FwTest)

# if __name__ == '__main__':
# 	t = [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23]
# 	Qo = [35.11, 47.77, 60.36, 76.29, 94.63, 103.44, 125.38, 147.39, 160.13, 145.11, 155.61, 137.93, 121.82, 109.27,
# 		  101.17, 98.40, 84.84, 83.86, 78.08, 69.18, 62.85, 58.89, 48.05]
# 	Np = [35.11, 82.88, 143.24, 219.53, 314.16, 417.60, 542.98, 690.38, 850.50, 995.62, 1151.22, 1289.15, 1410.97,
# 		  1520.24, 1621.41, 1719.82, 1804.66, 1888.52, 1966.59, 2035.77, 2098.62, 2157.51, 2205.56]
# 	Lp = [36.05, 85.27, 148.07, 227.98, 329.84, 446.89, 590.01, 765.93, 962.87, 1174.77, 1417.11, 1685.72, 1989.80,
# 		  2336.70, 2727.73, 3157.77, 3656.62, 4161.06, 4696.04, 5228.14, 5800.56, 6371.43, 6979.58]
# 	x = t[1:]
# 	y = Qo[1:]
# 	LpSample = Lp[6:]
# 	NpSample = Np[6:]
#
# 	alpha, beta, r2, b_fix = None, None, None, None
# 	if len(x) == len(y):
# 		beta, alpha, r2, b_fix = maxR2(x, y, len(x))
# 	else:
# 		print('length of x and y isn\'t equal.')
# 	print('alpha=', alpha, 'beta=', beta, 'r2=', r2, 'b=', b_fix)
# 	a = 10**alpha
# 	c = -1 / (2.303 * beta)
# 	Nr = a * c / (b_fix + 1)
# 	A, B, ABR2 = getAandB(NpSample, LpSample, len(NpSample))
# 	QoTest = getQoTest(t, len(t), a, b_fix, c)
# 	NpTest = getNpTest(t, len(t), a, b_fix, c)
# 	FwTest = getFwTest(Np, len(Np), A, B)
# 	Qmax, tm = getQmaxAndtm(a, b_fix, c)
#
